module fstats_smoothing use iso_fortran_env use fstats_errors use linalg, only : sort implicit none private public :: lowess contains ! ------------------------------------------------------------------------------ subroutine lowess(x, y, ys, fsmooth, nstps, del, rweights, resid) !! Computes the smoothing of a data set using a robust locally weighted !! scatterplot smoothing (LOWESS) algorithm. Fitted values are computed at !! each of the supplied x values. !! !! Remarks !! !! The code is a reimplementation of the LOWESS library. For a detailed !! understanding, see [this] !! (http://www.aliquote.org/cours/2012_biomed/biblio/Cleveland1979.pdf) !! paper by William Cleveland. real(real64), intent(in), dimension(:) :: x !! An N-element array containing the independent variable data. This !! array must be monotonically increasing. real(real64), intent(in), dimension(:) :: y !! An N-element array containing the dependent variable data. real(real64), intent(out), dimension(:) :: ys !! An N-element array where the smoothed results will be written. real(real64), intent(in), optional :: fsmooth !! An optional input that specifies the amount of smoothing. !! Specifically, this value is the fraction of points used to compute !! each value. As this value increases, the output becomes smoother. !! Choosing a value in the range of 0.2 to 0.8 typically results in a !! good fit. The default value is 0.2. integer(int32), intent(in), optional :: nstps !! An optional input that specifies the numb of iterations. If set to !! zero, a non-robust fit is returned. The default value is set to 2. real(real64), intent(in), optional :: del !! real(real64), intent(out), optional, dimension(:), target :: rweights !! An optional N-element array, that if supplied, will be used to !! return the weights given to each data point. real(real64), intent(out), optional, dimension(:), target :: resid !! An optional N-element array, that if supplied, will be used to !! return the residual. ! Parameters real(real64), parameter :: zero = 0.0d0 real(real64), parameter :: p2 = 2.0d-1 real(real64), parameter :: one = 1.0d0 real(real64), parameter :: three = 3.0d0 real(real64), parameter :: p001 = 1.0d-3 real(real64), parameter :: p999 = 0.999d0 ! Local Variables logical :: ok integer(int32) :: iter, i, j, nleft, nright, ns, last, m1, m2, n, nsteps real(real64) :: f, delta, d1, d2, denom, alpha, cut, eps, cmad, c1, c9, r real(real64), allocatable, target, dimension(:) :: rwdef, rsdef real(real64), pointer, dimension(:) :: rw, res ! Initialization n = size(x) if (present(fsmooth)) then f = fsmooth else f = p2 end if if (present(nstps)) then nsteps = nstps else nsteps = 2 end if if (present(del)) then delta = del else delta = 0.0d0 end if if (present(rweights)) then if (size(rweights) /= n) error stop FS_ARRAY_SIZE_ERROR rw => rweights else allocate(rwdef(n)) rw => rwdef end if if (present(resid)) then if (size(resid) /= n) error stop FS_ARRAY_SIZE_ERROR res => resid else allocate(rsdef(n)) res => rsdef end if ns = max(min(int(f * real(n), int32), n), 2) eps = epsilon(eps) ! Input Checking if (size(y) /= n) error stop FS_ARRAY_SIZE_ERROR if (size(ys) /= n) error stop FS_ARRAY_SIZE_ERROR ! Quick Return if (n < 2) then ys = y return end if ! Process do iter = 1, nsteps + 1 nleft = 1 nright = ns last = 0 i = 1 do do while (nright < n) d1 = x(i) - x(nleft) d2 = x(nright+1) - x(i) if (d1 <= d2) exit nleft = nleft + 1 nright = nright + 1 end do call lowest(x, y, x(i), ys(i), nleft, nright, res, iter > 1, & rw, ok) if (.not.ok) ys(i) = y(i) if (last < i - 1) then denom = x(i) - x(last) do j = last + 1, i - 1 alpha = (x(j) - x(last)) / denom ys(j) = alpha * ys(i) + (one - alpha) * ys(last) end do end if last = i cut = x(last) + delta do i = last + 1, n if (x(i) > cut) exit if (abs(x(i) - x(last)) < eps) then ys(i) = ys(last) last = i end if end do i = max(last + 1, i - 1) if (last >= n) exit end do res = y - ys if (iter > nsteps) exit rw = abs(res) call sort(rw, .true.) m1 = 1 + n / 2 m2 = n - m1 + 1 cmad = three * (rw(m1) + rw(m2)) c9 = p999 * cmad c1 = p001 * cmad do i = 1, n r = abs(res(i)) if (r <= c1) then rw(i) = one else if (r > c9) then rw(i) = zero else rw(i) = (one - (r / cmad)**2)**2 end if end do end do end subroutine ! ****************************************************************************** ! PRIVATE ROUTINES ! ------------------------------------------------------------------------------ ! REF: ! - https://en.wikipedia.org/wiki/Local_regression ! - http://www.aliquote.org/cours/2012_biomed/biblio/Cleveland1979.pdf subroutine lowest(x, y, xs, ys, nleft, nright, w, userw, rw, ok) ! Arguments real(real64), intent(in), dimension(:) :: x, y, rw ! N ELEMENT real(real64), intent(in) :: xs real(real64), intent(out) :: ys integer(int32), intent(in) :: nleft, nright real(real64), intent(out), dimension(:) :: w ! N ELEMENT logical, intent(in) :: userw logical, intent(out) :: ok ! Parameters real(real64), parameter :: zero = 0.0d0 real(real64), parameter :: one = 1.0d0 real(real64), parameter :: p001 = 1.0d-3 real(real64), parameter :: p999 = 0.999d0 ! Local Variables integer(int32) :: j, n, nrt real(real64) :: range, h, h9, h1, a, b, c, r ! Initialization n = size(x) range = x(n) - x(1) h = max(xs - x(nleft), x(nright) - xs) h9 = p999 * h h1 = p001 * h a = zero ! Process do j = nleft, n w(j) = zero r = abs(x(j) - xs) if (r <= h9) then if (r > h1) then w(j) = (one - (r / h)**3)**3 else w(j) = one end if if (userw) w(j) = rw(j) * w(j) a = a + w(j) else if (x(j) > xs) then exit end if end do nrt = j - 1 if (a <= zero) then ok = .false. else ok = .true. w(nleft:nrt) = w(nleft:nrt) / a if (h > zero) then a = zero do j = nleft, nrt a = a + w(j) * x(j) end do b = xs - a c = zero do j = nleft, nrt c = c + w(j) * (x(j) - a)**2 end do if (sqrt(c) > p001 * range) then b = b / c do j = nleft, nrt w(j) = w(j) * (one + b * (x(j) - a)) end do end if end if ys = zero do j = nleft, nrt ys = ys + w(j) * y(j) end do end if end subroutine ! ------------------------------------------------------------------------------ end module